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Summary. We present the phase diagram of a far from equilibrium system, mapped 
by Monte Carlo simulation. The model is a lattice gas of two species and holes. The 
two species are biased to hop in opposite directions and interact via excluded volume 
and nearest neighbor attractions. Three phases are found as function of temperature 
and charge density. 

Introduction. A comprehensive theory of nonequilibrium steady states is 
not yet available. This motivates the study of simple models for which analytic 
results can be compared to simulation data. While such models may seem far 
removed from physical phenomena, it is clearly easier to develop the tools for 
treating nonequilibrium systems in such a controlled setting before turning to 
the more demanding task of describing systems found in the natural world. 
For example, it is not yet known to what extent the idea of universality can be 
applied to nonequilibrium systems, or for that matter how to rigorously define 
a phase transition far from equilibrium. For the time being, we deal with these 
questions intuitively, amassing detailed descriptions of nonequilibrium systems 
while awaiting a general theory. In this paper we present the results of a Monte 
Carlo (MC) study of an interacting lattice gas driven far from equilibrium. 
In the first part we review the earlier work which motivates our study, and 
in the second part we describe the microscopic model and order parameters. 
Finally we describe the main results of our study, a phase diagram £P and a 
detailed finite-size scaling analysis of a particular limit of the model. 

Twenty years ago, Katz et al. introduced a simple modification (the 
"KLS model") of the attractive Ising lattice gas 3 , in which particle- hole 
exchanges along a particular lattice direction (which we call y) are coupled 
to a bias (E) which favors (suppresses) particle moves in the +y (—y) direc- 
tion. With the addition of periodic boundary conditions in the y direction, 
a nonzero particle current is established and the steady state distribution is 
not proportional to the Boltzmann distribution. A review of this work can be 
found in 4 . Simulations [7] and field-theoretic renormalization group 
studies 8 indicate that the model retains a continuous transition into a low- 
temperature phase separated state, though the universal behavior is distinct 
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from the Ising class (see below). Further, in contrast to the Ising model, the 
hopping bias allows only interfaces parallel to y. 

Now consider a generalization of such driven Ising models |U] to two species 
of particles which react in opposite senses to the bias. One type (+) is pushed 
in the +y direction while the other (— ) is pushed in the — y direction. In the 
limit E,T — > oo, while E/T = E finite the only interaction between particles 
is excluded volume. In this case a phase transition is found as the mass density 
(m) and E are varied. At low m and E the system is disordered. At sufficiently 
high m and E, the two species lock into a 'traffic jam': each blocks the other, 
and one observes a high-density strip perpendicular to y. 

If the high T, E constraint is lifted in the two-species model, then the 
attractive interactions will become significant over some range of T . This in- 
troduces the interesting possibility of transitions between the two types of 
order mentioned above. Imagine, e.g., beginning with the KLS model at half- 
filling. Now stay at half filling, but lower the charge density, i.e., change a 
few particles into ' — '. At some critical charge density, q c (E), the blocking 
transition may become stable. This possibility is investigated in the following 
sections by mapping the parameter space with MC simulations. Along the 
way we will find that it is quite difficult to find q c (E), as we lack a detailed 
understanding of the appropriate scaling forms in this region of parameter 
space. The subtleties of scaling arguments far from equilibrium will be illus- 
trated in the last section, where we will determine T c for the KLS model at 
finite E. 

Microscopic Model and Order Parameters. A configuration of the model 
is specified by a set of occupation variables, {s (r)}, where r = (x,y) labels a 
site on a fully periodic square lattice of dimensions L x x L y , and each s (r) 
can take the values +1, —1, or for a positive particle, negative particle, 
or hole. We also introduce the mass variable n (r) = |s (r)|. We will remain 
always at half-filling. The charge density is defined as q = L 1 L ^ r s (r). All 
particles interact via the usual Ising Hamiltonian, H = — 4J^ r r ' n ( r ) n ( r ')i 
independent of charge, with J > and the sum over nearest neighbors. J = 1 
is chosen arbitrarily: this merely sets an energy scale. The bias, E, points in 
the positive y-direction and is measured in units of J. A configuration evolves 
by selecting a nearest-neighbor bond at random; if occupied by a particle- 
hole pair, its contents are exchanged according to the Metropolis ^U] rate 
min {1, exp[— (AH — 5yEs(r))/T]}. The second term models the effect of the 
drive: if the particle, of charge s, is initially located at r, 5y is the change in its 
y-coordinate due to the jump. Thus, positive (negative) charges jump prefer- 
entially along (against) the field direction. The parameter T ("temperature") 
models the coupling to a thermal bath. The natural control parameters for 
our study are temperature T (measured in units of the Onsager value), the 
drive E and the charge density q. 

Conservation laws (for q and m) ensure spatially inhomogeneous ordered 
phases. We therefore select an order parameter sensitive to such structures, 
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Fig. 1. Phase diagram for E — 2 in charge density q and temperature T, in units 
of T c (E = 0). (+) denotes the first-order line, while the diamonds and x's denote 
the continuous lines. Configurations labelled: DO for disordered, HS for horizontal 
strip, VS for vertical strip 

i.e., the equal-time structure factor associated with the particle distribution 



(0(m x ,m y )) = ^ ^-^ r n(r)e' k - r \, k = 2tt(^,^). (•} denotes a con- 

figurational average, and the integers m x ,m y index the wave vector. For a 
perfect strip along the y-direction, {0(1, 0)} = 1 while a random configuration 
gives (0) = 0( L X L ). Except where noted, all simulations are run on 40 x 40 
lattices, starting from random initial configurations. One MC step (MCS) is 
defined as 2L x L y update attempts. The first 2 x 10 5 MCS are discarded, and 
measurements are taken every 200 MCS for the next 8 x 10 5 MCS. 

Monte Carlo Results: Phase Diagram. Although most detailed studies of 
the KLS model have considered the limit of infinite E to accentuate nonequi- 
librium properties, we are forced to consider finite values of E. As our two- 
species simulation stumbles through phase space it may find itself in a blocked 
configuration in a region of parameter space where such states should only 
be metastable. At large values of E this metastable configuration will persist 
far beyond the time of our simulation. Initial runs indicated that metastable 
lifetimes at E = 2 are reasonable. Note that at this value, jumps favoring E 
provide only one-half the change in energy required to break a single nearest- 
neighbor bond. It is therefore natural to wonder whether the transition at 
such small E remains in the KLS universality class. This question will be 
answered in the affirmative in the next two sections. 
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Fig. 1 shows the phase diagram in q and T for E = 2. The phase bound- 
aries are mapped by sweeping in T at fixed q while monitoring the order 
parameters (#(1, 0)) and (#(0, 1)) and their fluctuations. The far left point is 
the KLS transition at finite E. Moving right, parallel to the horizontal axis, 
we are changing + particles into — ; as long as the — are insufficient to form 
a blockade the transition from disorder into a parallel ("vertical") strip ap- 
parently remains continuous, as indicated by a peak in the fluctuations of 
(#(1,0)). We suspect that this transition can be described by the KLS field 
theory, with the additional complication of a small concentration of randomly 
distributed (but annealed) impurities. At the far right of the phase diagram 
the system contains equal numbers of + and — particles. There, we encounter 
a transition into a blocked phase ( "horizontal strip" ) which also appears con- 
tinuous, indicated here by a peak in the fluctuations of (#(0,1)). At larger 
values of q and lower T, we are able to observe a transition from the blocked 
phase into a vertical strip phase. This transition appears to be first-order, dis- 
playing hysteresis and metastability. The most interesting region in the phase 
diagram is where the three lines join at q c (E). Future studies will focus on 
scaling properties in the vicinity of this nonequilibrium bicritical point. 

Anisotropic Finite-Size Scaling. We now turn to a detailed discussion of 
the KLS transition at finite E. As mentioned in the introduction, a Langevin 
equation for a mesoscopic version of the local spin variable has been studied 
in great detail[H]- The most remarkable prediction of the field theory is a non- 
trivial anisotropy exponent, A > 1, so that wavevectors scale as kn ~ > 
physically this implies that domains of correlated spins grow faster in the 
field direction. In order to control finite size corrections, it is then necessary 
to account for this anisotropy The anisotropy introduces different corre- 
lation length exponents parallel(^n) and perpendicular {v±) to the drive, an 
effect which we will refer to as strong anisotropy, in contrast with, e.g., an 
Ising model with anisotropic interactions. The values of the critical expo- 
nents are known from an RG analysis to all orders; their values in d = 2 are 
(3 = 1/2,^|| = 3/2, v± = 1/2, A = 2|S]. Phenomenological scaling forms [7j 
involve two length scales, and L±, so that the order parameter scales as 

m(t,L\\,L±) = L|| ^""m (tL^ U{> , L"^" ]] and the scaling function de- 

pends on a "shape factor," S = L", x L^ 1 ^],^. Increasing the system size 
while holding S fixed allows us to approach T c without cutting off parallel 
correlations before transverse ones. We then use the (predicted) exponents to 
analyze our data. The validity of this approach will be judged by the quality 
of data collapse for m. In this way we will determine T c (E — 2), as it is the 
only fit parameter. Detailed work on carefully defined correlation functions 
and lengths [5] is forthcoming. 

Monte Carlo Results for E=2. We choose as our order parameter m = 
^ L n L J2 X y n{x, y^e 2 "^ 1 ^™-^! 1 ^^ \ since it is subject to smaller fluctuations 
than the structure factor. Most runs last for 1.2 x 10 6 MCS, though in larger 
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Fig. 2. Scaling of m at fixed S = .1575 for three system sizes: 24 x 54 (squares), 
28 x 86 (x) , 32 x 128 (triangles). The left (right) plot shows E = 20 [E = 2). The 
upper (lower) branch corresponds to T < T c (T > T c ). 

systems near criticality runs of 4.8 x 10 6 MCS were needed to ensure good 
statistics. The first 0.2 x 10 6 MCS were discarded and measurements were 
taken every 400 MCS thereafter. 

Fig. 2 shows the scaling of the order parameter for three different system 
sizes with S — .1575, for both E = 20 (effectively infinite) and E = 2. From 
these plots we estimate T c [E = 2) = 1.20(2). Notice the systematic deviations 
from scaling in the T < T c branch, which may be due to a small critical region 
or corrections to scaling from the marginal operator. Though we have not 
investigated this anomaly in detail, it occurs in the E = oo model as welljH] 
and has been observed in other nonequilibrium Ising models [TT). Elsewhere, 
the data collapse is of the same quality as in the E = oo case. We therefore 
have no reason to believe that the finite E transition falls into a different 
universality class. 

Conclusions. We have mapped out a slice of phase space for an interacting 
lattice gas of two species, driven far from equilibrium by a bias which drives 
a particle current. The phase diagram has two continuous lines which meet 
a first-order line at a critical charge density q c (E). In order to make more 
definitive claims, we need knowledge of scaling forms in the vicinity of the bi- 
critical point. It would be quite interesting to investigate this scaling behavior 
in order to learn about such points far from equilibrium. 
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